Machine Learning for Time Series Data in Python
Time series data is ubiquitous. Whether it be stock market fluctuations, sensor data recording climate change, or activity in the brain, any signal that changes over time can be described as a time series.
Machine learning has emerged as a powerful method for leveraging complexity in data in order to generate predictions and insights into the problem one is trying to solve. This course is an intersection between these two worlds of machine learning and time series data, and covers feature engineering, spectograms, and other advanced techniques in order to classify heartbeat sounds and predict stock prices.
PREREQUISITES:Manipulating Time Series Data in Python, Visualizing Time Series Data in Python, Machine Learning with scikit-learn
Download Datasets and Presentation slides for this post HERE
This post is 5th part of Time Series with Python Track from DataCamp. The track includes:
import numpy as np
import warnings
import pandas as pd
pd.set_option('display.expand_frame_repr', False)
import matplotlib.pyplot as plt
plt.style.use('seaborn')
import matplotlib as mpl
mpl.rcParams['figure.dpi'] = 150
mpl.rcParams['figure.figsize'] = (25, 20)
mpl.rcParams['axes.grid'] = True
warnings.filterwarnings("ignore")
Identifying a time series
Which of the following data sets is not considered time series data?
- Test grades for the last fall and spring semesters of high-school students.
- A student's attendance record each week of the semester.
- The school's national annual ranking since 2000.
- A list of the average length of each class at the school ---> You don't have timestamps for each data point, so it is not
Plotting a time series (I)
In this exercise, you'll practice plotting the values of two time series without the time component.
Two DataFrames, data and data2 are available in your workspace.
Instructions:
- Plot the values column of both the data sets on top of one another, one per axis object.
data = pd.read_csv('datasets/data.csv', index_col=[0])
data2 = pd.read_csv('datasets/data2.csv', index_col=[0])
display(data.head())
display(data2.head())
fig, axs = plt.subplots(2, 1, figsize=(15, 10))
data.iloc[:1000].plot(y='data_values', ax=axs[0], title="data")
data2.iloc[:1000].plot(y='data_values', ax=axs[1], title="data2")
plt.show()
Plotting a time series (II)
You'll now plot both the datasets again, but with the included time stamps for each (stored in the column called "time"). Let's see if this gives you some more context for understanding each time series data.
Instructions:
- Plot data and data2 on top of one another, one per axis object.
- The x-axis should represent the time stamps and the y-axis should represent the
data = pd.read_csv('datasets/data_with_tstamp.csv', index_col=[0])
data2 = pd.read_csv('datasets/data2_with_tstamp.csv', index_col=[0])
display(data.head())
display(data2.head())
fig, axs = plt.subplots(2, 1, figsize=(15, 10))
data.iloc[:1000].plot(x='time', y='data_values', ax=axs[0])
data2.iloc[:1000].plot(x='time', y='data_values', ax=axs[1])
plt.show()
Fitting a simple model:classification In this exercise, you'll use the iris dataset (representing petal characteristics of a number of flowers) to practice using the scikit-learn API to fit a classification model. You can see a sample plot of the data to the right.
Note: This course assumes some familiarity with Machine Learning and scikit-learn. For an introduction to scikit-learn, we recommend the Supervised Learning with Scikit-Learn and Preprocessing for Machine Learning in Python courses.
Instructions:
- Print the first five rows of data.
- Extract the "petal length (cm)" and "petal width (cm)" columns of data and assign it to X.
- Fit a model on X and y.
data = pd.read_csv('datasets/iris.csv', index_col=[0])
import matplotlib as mpl
mpl.rcParams.update(mpl.rcParamsDefault)
plt.style.use('default')
display(data.head())
print("data.target.unique: ",data.target.unique())
colors = np.where(data["target"]==1,'purple','red')
data.plot.scatter(x='petal length (cm)', y='petal width (cm)', c=colors)
from sklearn.svm import LinearSVC
# Construct data for the model
X = data[['petal length (cm)','petal width (cm)']]
y = data[['target']]
# Fit the model
model = LinearSVC()
model.fit(X, y)
Predicting using a classification model
Now that you have fit your classifier, let's use it to predict the type of flower (or class) for some newly-collected flowers.
Information about petal width and length for several new flowers is stored in the variable targets. Using the classifier you fit, you'll predict the type of each flower.
Instructions:
- Predict the flower type using the array X_predict.
- Run the given code to visualize the predictions.
targets = pd.read_csv('datasets/targets.csv', index_col=[0])
X_predict = targets[['petal length (cm)', 'petal width (cm)']]
# Predict with the model
predictions = model.predict(X_predict)
print(predictions)
# [2 2 2 1 1 2 2 2 2 1 2 1 1 2 1 1 2 1 2 2]
# Visualize predictions and actual values
plt.scatter(X_predict['petal length (cm)'], X_predict['petal width (cm)'],
c=predictions, cmap=plt.cm.coolwarm)
plt.title("Predicted class values")
plt.show()
""" Note that the output of your predictions are all integers, representing that datapoint's predicted class."""
Fitting a simple model:regression
In this exercise, you'll practice fitting a regression model using data from the California housing market. A DataFrame called housing is available in your workspace. It contains many variables of data (stored as columns). Can you find a relationship between the following two variables?
- "MedHouseVal": the median house value for California districts (in $100,000s of dollars)
- "AveRooms" : average number of rooms per dwelling
Instructions:
- Prepare X and y DataFrames using the data in housing.
- X should be the Median House Value, y average number of rooms per dwelling.
- Fit a regression model that uses these variables (remember to shape the variables correctly!).
- Don't forget that each variable must be the correct shape for scikit-learn to use it!
housing = pd.read_csv("datasets/housing.csv", index_col=[0])
display(housing.info())
housing.plot.scatter(y='AveRooms', x='MedHouseVal')
housing[housing.isna().any(axis=1)]
housing.dropna(inplace=True)
from sklearn import linear_model
# Prepare input and output DataFrames
X = housing[['MedHouseVal']]
y = housing[['AveRooms']]
# Fit the model
model = linear_model.LinearRegression()
model.fit(X,y)
"""In regression, the output of your model is a continuous array of numbers, not class identity."""
Predicting using a regression model
Now that you've fit a model with the California housing data, lets see what predictions it generates on some new data. You can investigate the underlying relationship that the model has found between inputs and outputs by feeding in a range of numbers as inputs and seeing what the model predicts for each input.
A 1-D array new_inputs consisting of 100 "new" values for "MedHouseVal" (median house value) is available in your workspace along with the model you fit in the previous exercise.
Instructions:
- Review new_inputs in the shell.
- Reshape new_inputs appropriately to generate predictions.
- Run the given code to visualize the predictions.
new_inputs_df = pd.read_csv('./datasets/newinputs.csv', index_col=[0])
new_inputs = new_inputs_df.to_numpy()
new_inputs.reshape((100,))
predictions = model.predict(new_inputs.reshape(-1,1))
# Visualize the inputs and predicted values
plt.scatter(new_inputs, predictions, color='r', s=3)
plt.xlabel('inputs')
plt.ylabel('predictions')
plt.show()
""" Here the red line shows the relationship that your model found. As the number of rooms grows, the median house value rises linearly."""
Inspecting the classification data
In these final exercises of this chapter, you'll explore the two datasets you'll use in this course.
The first is a collection of heartbeat sounds. Hearts normally have a predictable sound pattern as they beat, but some disorders can cause the heart to beat abnormally. This dataset contains a training set with labels for each type of heartbeat, and a testing set with no labels. You'll use the testing set to validate your models.
As you have labeled data, this dataset is ideal for classification. In fact, it was originally offered as a part of a public Kaggle competition.
Instructions:
- Use glob to return a list of the .wav files in data_dir directory.
- Import the first audio file in the list using librosa.
- Generate a time array for the data.
- Plot the waveform for this file, along with the time array.
import librosa as lr
from glob import glob
data_dir = './files'
# List all the wav files in the folder
audio_files = glob(data_dir + '/*.wav')
print(len(audio_files))
# Read in the first audio file, create the time array
audio, sfreq = lr.load(audio_files[0])
time = np.arange(0, len(audio)) / sfreq
# Plot audio over time
fig, ax = plt.subplots()
ax.plot(time, audio)
ax.set(xlabel='Time (s)', ylabel='Sound Amplitude')
plt.show()
Output:
Inspecting the regression data
The next dataset contains information about company market value over several years of time. This is one of the most popular kind of time series data used for regression. If you can model the value of a company as it changes over time, you can make predictions about where that company will be in the future. This dataset was also originally provided as part of a public Kaggle competition.
In this exercise, you'll plot the time series for a number of companies to get an understanding of how they are (or aren't) related to one another.
Instructions:
- Import the data with Pandas (stored in the file 'prices.csv').
- Convert the index of data to datetime.
- Loop through each column of data and plot the the column's values over time.
mpl.rcParams['figure.figsize'] = (15, 10)
data = pd.read_csv('./datasets/prices_new.csv', index_col=0)
# Convert the index of the DataFrame to datetime
data.index = pd.to_datetime(data.index)
print(data.head())
# Loop through each column, plot its values over time
fig, ax = plt.subplots()
for column in data.columns:
data[column].plot(ax=ax, label=column)
ax.legend()
plt.show()
""" Note that each company's value is sometimes correlated with others, and sometimes not. Also note there are a lot of 'jumps' in there - what effect do you think these jumps would have on a predictive model? """
Some preprocessing steps to load the data for next sections.
import csv
import os
import re
def convert_txt_to_csv(file_path):
new_file_name = os.path.splitext(file_path)[0] + '.csv'
with open(file_path, "r") as f:
plain_str = f.read()
print(plain_str[:10], plain_str[-10:])
csv_data = re.match(r".*?b\'(.*)\'.*?", plain_str).group(1) # remove byte traces
print(csv_data[:10], csv_data[-10:])
csv_data = csv_data.replace('\\n', '\n') # remove \\n with actual new line chars
with open(new_file_name, 'w') as out:
out.write(csv_data)
files_to_convert = [
'./datasets/abnormal.txt',
'./datasets/normal.txt'
]
for file_path in files_to_convert:
convert_txt_to_csv(file_path)
abnormal_df = pd.read_csv('./datasets/abnormal.csv', index_col=[0])
display(abnormal_df.tail())
print(abnormal_df.info())
# Read in the data
normal_df = pd.read_csv('./datasets/normal.csv', index_col=[0])
display(normal_df.tail())
print(normal_df.info())
normal = normal_df.to_numpy()
print(normal.shape)
abnormal = abnormal_df.to_numpy()
print(abnormal.shape)
sfreq = 2205
# import inspect
# source = inspect.getsource(show_plot_and_make_titles)
# print(source)
def show_plot_and_make_titles():
axs[0, 0].set(title="Normal Heartbeats")
axs[0, 1].set(title="Abnormal Heartbeats")
plt.tight_layout()
plt.show()
Many repetitions of sounds
In this exercise, you'll start with perhaps the simplest classification technique:averaging across dimensions of a dataset and visually inspecting the result. You'll use the heartbeat data described in the last chapter. Some recordings are normal heartbeat activity, while others are abnormal activity. Let's see if you can spot the difference.
Two DataFrames, normal and abnormal, each with the shape of (n_times_points, n_audio_files) containing the audio for several heartbeats are available in your workspace. Also, the sampling frequency is loaded into a variable called sfreq. A convenience plotting function show_plot_and_make_titles() is also available in your workspace.
Instructions:
- First, create the time array for these audio files (all audios are the same length).
- Then, stack the values of the two DataFrames together (
normalandabnormal, in that order) so that you have a single array of shape (n_audio_files,n_times_points). - Finally, use the code provided to loop through each list item / axis, and plot the audio over time in the corresponding axis object.
- You'll plot normal heartbeats in the left column, and abnormal ones in the right column
import matplotlib as mpl
mpl.rcParams['figure.dpi'] = 200
mpl.rcParams['figure.figsize'] = (30, 15)
fig, axs = plt.subplots(3, 2, figsize=(30, 15), sharex=True, sharey=True)
# Calculate the time array
time = np.arange(normal.shape[0]) / sfreq
# Stack the normal/abnormal audio so you can loop and plot
stacked_audio = np.hstack([normal, abnormal]).T
# Loop through each audio file / ax object and plot
# .T.ravel() transposes the array, then unravels it into a 1-D vector for looping
for iaudio, ax in zip(stacked_audio, axs.T.ravel()):
ax.plot(time, iaudio)
show_plot_and_make_titles()
""" As you can see there is a lot of variability in the raw data, let's see if you can average out some of that noise to notice a difference."""
Invariance in time
While you should always start by visualizing your raw data, this is often uninformative when it comes to discriminating between two classes of data points. Data is usually noisy or exhibits complex patterns that aren't discoverable by the naked eye.
Another common technique to find simple differences between two sets of data is to average across multiple instances of the same class. This may remove noise and reveal underlying patterns (or, it may not).
In this exercise, you'll average across many instances of each class of heartbeat sound.
The two DataFrames (normal and abnormal) and the time array (time) from the previous exercise are available in your workspace.
Instructions:
- Average across the audio files contained in normal and abnormal, leaving the time dimension.
- Visualize these averages over time.
files_to_convert = [
'./datasets/invariance_abnormal.txt',
'./datasets/invariance_normal.txt'
]
for file_path in files_to_convert:
convert_txt_to_csv(file_path)
normal = pd.read_csv('./datasets/invariance_normal.csv', index_col=[0])
abnormal = pd.read_csv('./datasets/invariance_abnormal.csv', index_col=[0])
mean_normal = np.mean(normal, axis=1)
mean_abnormal = np.mean(abnormal, axis=1)
# Plot each average over time
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(17, 5), sharey=True)
ax1.plot(time, mean_normal)
ax1.set(title="Normal Data")
ax2.plot(time, mean_abnormal)
ax2.set(title="Abnormal Data")
plt.show()
"""Do you see a noticeable difference between the two? Maybe, but it's quite noisy. Let's see how you can dig into the data a bit further."""
Build a classification model
While eye-balling differences is a useful way to gain an intuition for the data, let's see if you can operationalize things with a model. In this exercise, you will use each repetition as a datapoint, and each moment in time as a feature to fit a classifier that attempts to predict abnormal vs. normal heartbeats using only the raw data.
We've split the two DataFrames (normal and abnormal) into X_train, X_test, y_train, and y_test.
Instructions:
- Create an instance of the Linear SVC model and fit the model using the training data.
- Use the testing data to generate predictions with the model.
- Score the model using the provided code.
files_to_convert = [ './datasets/abnormal_full.txt',
'./datasets/normal_full.txt' ]
for file_path in files_to_convert:
convert_txt_to_csv(file_path)
normal = pd.read_csv('./datasets/normal_full.csv', index_col=0).select_dtypes(np.float64).astype(np.float32)
abnormal = pd.read_csv('./datasets/abnormal_full.csv', index_col=0).select_dtypes(np.float64).astype(np.float32)
display(normal.head())
display(abnormal.head())
print(normal.shape)
print(abnormal.shape)
stacked_audio = np.hstack([normal, abnormal]).T
X = stacked_audio
X.shape
y = np.concatenate((np.full((len(normal.T), 1), 'normal'), np.full((len(abnormal.T), 1), 'abnormal')))
y.shape
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42)
print(type(X), X[:2],X.shape)
print(type(y),y[:2], y.shape)
len(y_test)
from sklearn.svm import LinearSVC
# Initialize and fit the model
model = LinearSVC()
model.fit(X_train,y_train)
# Generate predictions and score them manually
predictions = model.predict(X_test)
print(sum(predictions == y_test.squeeze()) / len(y_test))
"""Note that your predictions didn't do so well. That's because the features you're using as inputs to the model (raw data) aren't very good at differentiating classes. Next, you'll explore how to calculate some more complex features that may improve the results."""
Calculating the envelope of sound
One of the ways you can improve the features available to your model is to remove some of the noise present in the data. In audio data, a common way to do this is to smooth the data and then rectify it so that the total amount of sound energy over time is more distinguishable. You'll do this in the current exercise.
A heartbeat file is available in the variable audio.
Instructions:
- Visualize the raw audio you'll use to calculate the envelope.
import h5py
path = './datasets/audio_munged.hdf5'
f = h5py.File(path,'r')
list(f.keys())
f['h5io'].keys()
sfreq = f['h5io']['key_sfreq'][0] # 2205
key_data = pd.read_hdf(path, key='h5io/key_data')
audio = key_data.squeeze()
display(audio.shape)
key_meta = pd.read_hdf(path, key='h5io/key_meta')
display(key_meta.head())
display(key_meta.info())
labels = key_meta["label"].to_numpy()
labels
Opps! The dataset that Datacamp provided could not be loaded like the one in their exercise. Let's load it manually.
audio = pd.read_csv('./datasets/audio.csv', index_col=0)
audio.plot(figsize=(10, 5))
plt.show()
- Rectify the audio.
- Plot the result.
audio_rectified = audio.apply(np.abs)
# Plot the result
audio_rectified.plot(figsize=(10, 5))
plt.show()
- Smooth the audio file by applying a rolling mean.
- Plot the result.
audio_rectified_smooth = audio_rectified.rolling(50).mean()
# Plot the result
audio_rectified_smooth.plot(figsize=(10, 5))
plt.show()
By calculating the envelope of each sound and smoothing it, you've eliminated much of the noise and have a cleaner signal to tell you when a heartbeat is happening.
Calculating features from the envelope
Now that you've removed some of the noisier fluctuations in the audio, let's see if this improves your ability to classify.
audio_rectified_smooth from the previous exercise is available in your workspace.
Instructions:
- Calculate the mean, standard deviation, and maximum value for each heartbeat sound.
- Column stack these stats in the same order.
- Use cross-validation to fit a model on each CV iteration.
path = './datasets/audio_munged.hdf5'
df = pd.read_hdf(path, key='h5io/key_data')
audio = df.squeeze()
# Smooth by applying a rolling mean
audio_rectified = audio.apply(np.abs)
audio_rectified_smooth = audio_rectified.rolling(50).mean()
audio_rectified_smooth
model = LinearSVC()
means = np.mean(audio_rectified_smooth, axis=0)
stds = np.std(audio_rectified_smooth, axis=0)
maxs = np.max(audio_rectified_smooth, axis=0)
# Create the X and y arrays
X = np.column_stack([means, stds, maxs])
y = labels.reshape([-1, 1])
# Fit the model and score on testing data
from sklearn.model_selection import cross_val_score
percent_score = cross_val_score(model, X, y, cv=5)
print(np.mean(percent_score))
This model is both simpler (only 3 features) and more understandable (features are simple summary statistics of the data).
Derivative features:The tempogram One benefit of cleaning up your data is that it lets you compute more sophisticated features. For example, the envelope calculation you performed is a common technique in computing tempo and rhythm features. In this exercise, you'll use
librosato compute some tempo and rhythm features for heartbeat data, and fit a model once more.
Note that librosa functions tend to only operate on numpy arrays instead of DataFrames, so we'll access our Pandas data as a Numpy array with the .values attribute.
Instructions:
- Use
librosato calculate a tempogram of each heartbeat audio. - Calculate the mean, standard deviation, and maximum of each tempogram (this time using DataFrame methods)
- Column stack these tempo features (mean, standard deviation, and maximum) in the same order.
- Score the classifier with cross-validation.
model = LinearSVC()
import librosa as lr
# Calculate the tempo of the sounds
tempos = []
for col, i_audio in audio.items():
tempos.append(lr.beat.tempo(i_audio.values, sr=sfreq, hop_length=2**6, aggregate=None))
# Convert the list to an array so you can manipulate it more easily
tempos = np.array(tempos)
# Calculate statistics of each tempo
tempos_mean = tempos.mean(axis=-1)
tempos_std = tempos.std(axis=-1)
tempos_max = tempos.max(axis=-1)
X = np.column_stack([means, stds, maxs, tempos_mean, tempos_std, tempos_max])
y = labels.reshape([-1, 1])
# Fit the model and score on testing data
percent_score = cross_val_score(model, X, y, cv=5)
print(np.mean(percent_score))
"""Note that your predictive power may not have gone up (because this dataset is quite small), but you now have a more rich feature representation of audio that your model can use!"""
Spectrograms of heartbeat audio
Spectral engineering is one of the most common techniques in machine learning for time series data. The first step in this process is to calculate a spectrogram of sound. This describes what spectral content (e.g., low and high pitches) are present in the sound over time. In this exercise, you'll calculate a spectrogram of a heartbeat audio file.
We've loaded a single heartbeat sound in the variable audio.
Instructions:
- Import the short-time fourier transform (
stft) function fromlibrosa.core. - Calculate the spectral content (using the short-time fourier transform function) of
audio. - Convert the spectogram (
spec) to decibels. - Visualize the spectogram.
key_data = pd.read_hdf(path, key='h5io/key_data')
audio = key_data.iloc[:, 0].to_numpy()
print(type(audio), audio.shape)
audio
time = key_data.index.to_numpy()
time.shape
from librosa.core import stft
# Prepare the STFT
HOP_LENGTH = 2**4
SIZE_WINDOW = 2**7
spec = stft(audio, hop_length=HOP_LENGTH, n_fft= SIZE_WINDOW)
spec
times_spec = time[::HOP_LENGTH]
from librosa.core import amplitude_to_db
from librosa.display import specshow
# Convert into decibels
spec_db = amplitude_to_db(spec)
# Compare the raw audio to the spectrogram of the audio
fig, axs = plt.subplots(2, 1, figsize=(15, 15), sharex=True)
axs[0].plot(time, audio)
specshow(spec_db, sr=sfreq, x_axis='time', y_axis='hz', hop_length=HOP_LENGTH)
plt.show()
""" Do you notice that the heartbeats come in pairs, as seen by the vertical lines in the spectrogram? """
Engineering spectral features
As you can probably tell, there is a lot more information in a spectrogram compared to a raw audio file. By computing the spectral features, you have a much better idea of what's going on. As such, there are all kinds of spectral features that you can compute using the spectrogram as a base. In this exercise, you'll look at a few of these features.
The spectogram spec from the previous exercise is available in your workspace.
Instructions:
- Calculate the spectral bandwidth as well as the spectral centroid of the spectrogram by using functions in
librosa.feature.
spec = np.absolute(spec.astype(np.float32))
spec
import librosa as lr
# Calculate the spectral centroid and bandwidth for the spectrogram
bandwidths = lr.feature.spectral_bandwidth(S=spec)[0]
centroids = lr.feature.spectral_centroid(S=spec)[0]
times_spec = time[::HOP_LENGTH]
import matplotlib as mpl
mpl.rcParams['figure.dpi'] = 300
from librosa.core import amplitude_to_db
from librosa.display import specshow
# Convert spectrogram to decibels for visualization
spec_db = amplitude_to_db(spec)
# Display these features on top of the spectrogram
fig, ax = plt.subplots(figsize=(15, 10))
specshow(spec_db, ax = ax, x_axis='time', y_axis='hz', hop_length=HOP_LENGTH)
ax.plot(times_spec, centroids)
ax.fill_between(times_spec, centroids - bandwidths / 2, centroids + bandwidths / 2, alpha=.5)
ax.set(ylim=[None, 6000])
plt.show()
"""As you can see, the spectral centroid and bandwidth characterize the spectral content in each sound over time. They give us a summary of the spectral content that we can use in a classifier."""
Combining many features in a classifier
You've spent this lesson engineering many features from the audio data - some contain information about how the audio changes in time, others contain information about the spectral content that is present.
The beauty of machine learning is that it can handle all of these features at the same time. If there is different information present in each feature, it should improve the classifier's ability to distinguish the types of audio. Note that this often requires more advanced techniques such as regularization, which we'll cover in the next chapter.
For the final exercise in the chapter, we've loaded many of the features that you calculated before. Combine all of them into an array that can be fed into the classifier, and see how it does.
Instructions:
- Loop through each spectrogram, calculating the mean spectral bandwidth and centroid of each.
from librosa.core import stft
import librosa as lr
# Prepare the STFT
HOP_LENGTH = 2**4
SIZE_WINDOW = 2**7
path = './datasets/audio_munged.hdf5'
key_data = pd.read_hdf(path, key='h5io/key_data')
cols_list = [col_values.to_numpy() for _, col_values in key_data.iteritems()]
def mySTFT(col_data, H_L =HOP_LENGTH, S_Z = SIZE_WINDOW):
col_spec = stft(col_data, hop_length= H_L, n_fft= S_Z)
return np.abs(col_spec.astype(np.float32))
# mySTFT(cols_list[0])
spectrograms = [ mySTFT(col) for col in cols_list]
spectrograms[0].shape
bandwidths = []
centroids = []
for spec in spectrograms:
# Calculate the mean spectral bandwidth
this_mean_bandwidth = np.mean(lr.feature.spectral_bandwidth(S=spec))
# Calculate the mean spectral centroid
this_mean_centroid = np.mean(lr.feature.spectral_centroid(S=spec))
# Collect the values
bandwidths.append(this_mean_bandwidth)
centroids.append(this_mean_centroid)
path = './datasets/audio_munged.hdf5'
df = pd.read_hdf(path, key='h5io/key_data')
audio = df.squeeze()
key_meta = pd.read_hdf(path, key='h5io/key_meta')
labels = key_meta["label"].to_numpy()
# Smooth by applying a rolling mean
audio_rectified = audio.apply(np.abs)
audio_rectified_smooth = audio_rectified.rolling(50).mean()
# Calculate stats
means = np.mean(audio_rectified_smooth, axis=0)
stds = np.std(audio_rectified_smooth, axis=0)
maxs = np.max(audio_rectified_smooth, axis=0)
# Calculate the tempo of the sounds
tempos = []
for col, i_audio in audio.items():
tempos.append(lr.beat.tempo(i_audio.values, sr=sfreq, hop_length=2**6, aggregate=None))
# Convert the list to an array so you can manipulate it more easily
tempos = np.array(tempos)
# Calculate statistics of each tempo
tempos_mean = tempos.mean(axis=-1)
tempos_std = tempos.std(axis=-1)
tempos_max = tempos.max(axis=-1)
from sklearn.svm import LinearSVC
from sklearn.model_selection import cross_val_score
model = LinearSVC()
- Column stack all the features to create the array X.
- Score the classifier with cross-validation.
X = np.column_stack([means, stds, maxs, tempos_mean, tempos_max, tempos_std, bandwidths, centroids])
y = labels.reshape([-1, 1])
# Fit the model and score on testing data
percent_score = cross_val_score(model, X, y, cv=5)
print(np.mean(percent_score))
Oh, the result does not looks like the DataCamp's output. I checked, and it looks like the labels that they provide does not match with what was using. Here is the label raw data.
labels = np.array(['normal', 'normal', 'normal', 'murmur', 'normal', 'normal',
'normal', 'murmur', 'normal', 'murmur', 'normal', 'normal',
'normal', 'murmur', 'murmur', 'normal', 'normal', 'murmur',
'murmur', 'normal', 'murmur', 'murmur', 'murmur', 'murmur',
'normal', 'normal', 'murmur', 'normal', 'normal', 'murmur',
'murmur', 'murmur', 'murmur', 'murmur', 'murmur', 'normal',
'normal', 'murmur', 'murmur', 'murmur', 'normal', 'murmur',
'murmur', 'normal', 'normal', 'normal', 'murmur', 'murmur',
'murmur', 'normal', 'normal', 'normal', 'normal', 'murmur',
'normal', 'normal', 'murmur', 'murmur', 'murmur', 'murmur'],
dtype=object)
Now we fit the model and score it again.
model3 = LinearSVC()
# Create X and y arrays
X = np.column_stack([means, stds, maxs, tempos_mean, tempos_max, tempos_std, bandwidths, centroids])
y = labels.reshape([-1, 1])
# Fit the model and score on testing data
percent_score = cross_val_score(model3, X, y, cv=5)
print(np.mean(percent_score))
""" There we go, this was the shown result. Note that if we rerun the cell, the score's result will change. """
You calculated many different features of the audio, and combined each of them under the assumption that they provide independent information that can be used in classification. You may have noticed that the accuracy of your models varied a lot when using different set of features. This chapter was focused on creating new "features" from raw data and not obtaining the best accuracy. To improve the accuracy, you want to find the right features that provide relevant information and also build models on much larger data.
Predicting Time Series Data
If you want to predict patterns from data over time, there are special considerations to take in how you choose and construct your model. This chapter covers how to gain insights into the data before fitting your model, as well as best-practices in using predictive modeling for time series data.
Introducing the dataset
As mentioned in the video, you'll deal with stock market prices that fluctuate over time. In this exercise you've got historical prices from two tech companies (Ebay and Yahoo) in the DataFrame prices. You'll visualize the raw data for the two companies, then generate a scatter plot showing how the values for each company compare with one another. Finally, you'll add in a "time" dimension to your scatter plot so you can see how this relationship changes over time.
The data has been loaded into a DataFrame called prices.
Instructions:
- Plot the data in prices. Pay attention to any irregularities you notice.
- Generate a scatter plot with the values of Ebay on the x-axis, and Yahoo on the y-axis. Look up the symbols for both companies from the column names of the DataFrame.
- Finally, encode time as the color of each datapoint in order to visualize how the relationship between these two variables changes.
prices = pd.read_csv('./datasets/prices_ebay_yahoo.csv', index_col=0, parse_dates=True)
prices.info()
import matplotlib.pyplot as plt
plt.style.use('ggplot')
import matplotlib as mpl
mpl.rcParams['figure.dpi'] = 150
mpl.rcParams['figure.figsize'] = (15, 12)
prices.plot()
# add legend and set position to upper left
plt.legend( loc='upper left')
plt.show()
prices.plot.scatter('EBAY', 'YHOO')
plt.legend( loc='upper left')
plt.show()
prices.plot.scatter('EBAY', 'YHOO', c=prices.index, cmap=plt.cm.viridis, colorbar=False, sharex=False, fontsize=18, title='Color = time')
# Define your mappable for colorbar creation
sm = plt.cm.ScalarMappable(cmap='viridis',
norm=plt.Normalize(vmin=prices.index.min().value,
vmax=prices.index.max().value))
cbar = plt.colorbar(sm)
# Change the numeric ticks into ones that match the x-axis
cbar.ax.set_yticklabels(pd.to_datetime(cbar.get_ticks()).strftime(date_format='%b %Y'))
Let's try with a different view.
df = prices
fig,ax=plt.subplots()
sca = ax.scatter(df.index,df['EBAY'],c=df.index,cmap='viridis', s=10)
sca2 = ax.scatter(df.index,df['YHOO'],c=df.index,cmap='jet_r', s=10)
cbar = plt.colorbar(sca)
cbar.ax.set_title("EBAY")
cbar.ax.set_yticklabels(pd.to_datetime(cbar.get_ticks()).strftime(date_format='%b %Y'))
cbar2 = plt.colorbar(sca2)
cbar2.ax.set_title("YHOO")
cbar2.ax.set_yticklabels(pd.to_datetime(cbar.get_ticks()).strftime(date_format='%b %Y'))
fig.autofmt_xdate()
ax.set_title("Color = Time")
plt.show()
"""As you can see, these two time series seem somewhat related to each other, though its a complex relationship that changes over time."""
Fitting a simple regression model
Now we'll look at a larger number of companies. Recall that we have historical price values for many companies. Let's use data from several companies to predict the value of a test company. You'll attempt to predict the value of the Apple stock price using the values of NVidia, Ebay, and Yahoo. Each of these is stored as a column in the all_prices DataFrame. Below is a mapping from company name to column name:ebay: "EBAY"
nvidia: "NVDA"
yahoo: "YHOO"
apple: "AAPL"
We'll use these columns to define the input/output arrays in our model.
Instructions:
- Create the X and y arrays by using the column names provided.
- The input values should be from the companies "ebay", "nvidia", and "yahoo"
- The output values should be from the company "apple"
- Use the data to train and score the model with cross-validation.
import re, csv, os
file_path = './datasets/all_prices.txt'
with open(file_path, "r") as f:
plain_str = f.read()
print(plain_str[:10], plain_str[-10:])
csv_data = re.match(r".*?b\'(.*)\'.*?", plain_str).group(1) # remove byte traces
print("START:\n",csv_data[:250],"\nEND:", csv_data[-250:])
csv_data = csv_data.replace('\\n', '\n') # remove \\n with actual new line chars
print("START:\n",csv_data[:250],"\nEND:", csv_data[-250:])
with open(os.path.splitext(file_path)[0] + '.csv', 'w') as out:
out.write(csv_data)
all_prices = pd.read_csv('./datasets/all_prices.csv', index_col=0, parse_dates=True)
all_prices.head()
from sklearn.linear_model import Ridge
from sklearn.model_selection import cross_val_score
# Use stock symbols to extract training data
X = all_prices[['EBAY', 'NVDA', 'YHOO']]
y = all_prices[['AAPL']]
# Fit and score the model with cross-validation
scores = cross_val_score(Ridge(), X, y, cv=3)
print(scores)
""" As you can see, fitting a model with raw data doesn't give great results."""
Visualizing predicted values
When dealing with time series data, it's useful to visualize model predictions on top of the "actual" values that are used to test the model.
In this exercise, after splitting the data (stored in the variables X and y) into training and test sets, you'll build a model and then visualize the model's predictions on top of the testing data in order to estimate the model's performance.
Instructions:
- Split the data (X and y) into training and test sets.
- Use the training data to train the regression model.
- Then use the testing data to generate predictions for the model.
X = pd.read_csv('./datasets/X.csv', index_col=[0]).to_numpy()
y = pd.read_csv('./datasets/y.csv', index_col=[0]).to_numpy()
display(X.shape)
display(y.shape)
from sklearn.model_selection import train_test_split
from sklearn.metrics import r2_score
# Split our data into training and test sets
X_train, X_test, y_train, y_test = train_test_split(X, y,
train_size=.8, shuffle=False, random_state=1)
# Fit our model and generate predictions
model = Ridge()
model.fit(X_train, y_train)
predictions = model.predict(X_test)
score = r2_score(y_test, predictions)
print(score)
- Plot a time series of the predicted and "actual" values of the testing data.
mpl.rcParams['figure.figsize'] = (30, 20)
mpl.rcParams["xtick.labelsize"], mpl.rcParams["ytick.labelsize"] = 25,25
fig, ax = plt.subplots()
ax.plot(y_test, color='k', lw=3)
ax.plot(predictions, color='r', lw=2)
plt.show()
Visualizing messy data
Let's take a look at a new dataset - this one is a bit less-clean than what you've seen before.
As always, you'll first start by visualizing the raw data. Take a close look and try to find datapoints that could be problematic for fitting models.
The data has been loaded into a DataFrame called prices.
Instructions:
- Visualize the time series data using Pandas.
- Calculate the number of missing values in each time series. Note any irregularities that you can see. What do you think they are?
prices = pd.read_csv('./datasets/prices_ebay_yahoo_nvidia.csv', parse_dates=True, index_col=0)
prices
prices.plot(legend=True)
plt.tight_layout()
plt.show()
# Count the missing values of each time series
missing_values = prices.isna().sum()
print(missing_values)
""" In the plot, you can see there are clearly missing chunks of time in your data. There also seem to be a few 'jumps' in the data. How can you deal with this?"""
Imputing missing values
When you have missing data points, how can you fill them in?
In this exercise, you'll practice using different interpolation methods to fill in some missing values, visualizing the result each time. But first, you will create the function (interpolate_and_plot()) you'll use to interpolate missing data points and plot them.
A single time series has been loaded into a DataFrame called prices.
Instructions:
- Create a boolean mask for missing values and interpolate the missing values using the
interpolationargument of the function. - Interpolate using the latest non-missing value and plot the results. Recall that
interpolate_and_plot's second input is a string specifying the kind of interpolation to use. - Interpolate linearly and plot the results.
- Interpolate with a quadratic function and plot the results.
def interpolate_and_plot(prices, interpolation):
# Create a boolean mask for missing values
missing_values = prices.isna()
# Interpolate the missing values
prices_interp = prices.interpolate(interpolation)
# Plot the results, highlighting the interpolated values in black
fig, ax = plt.subplots(figsize=(10, 5))
prices_interp.plot(color='k', alpha=.6, ax=ax, legend=False)
# Now plot the interpolated values on top in red
prices_interp[missing_values].plot(ax=ax, color='r', lw=3, legend=False)
plt.show()
plt.style.use('seaborn')
mpl.rcParams['figure.figsize'] = (30, 20)
mpl.rcParams["xtick.labelsize"], mpl.rcParams["ytick.labelsize"] = 12,12
prices = prices.head(1278)
prices
interpolation_type = 'zero'
interpolate_and_plot(prices, interpolation_type)
interpolation_type = 'quadratic'
interpolate_and_plot(prices, interpolation_type)
""" Correct! When you interpolate, the pre-existing data is used to infer the values of missing data. As you can see, the method you use for this has a big effect on the outcome."""
Transforming raw data
In the last chapter, you calculated the rolling mean. In this exercise, you will define a function that calculates the percent change of the latest data point from the mean of a window of previous data points. This function will help you calculate the percent change over a rolling window.
This is a more stable kind of time series that is often useful in machine learning.
Instructions:
- Define a percent_change function that takes an input time series and does the following:
- Extract all but the last value of the input series (assigned to previous_values) and the only the last value of the timeseries ( assigned to last_value)
- Calculate the percentage difference between the last value and the mean of earlier values.
- Using a rolling window of 20, apply this function to prices, and visualize it using the given code.
prices = all_prices[["EBAY" , "NVDA" , "YHOO" , "AAPL"]]
prices
mpl.rcParams['figure.figsize'] = (30, 30)
mpl.rcParams["xtick.labelsize"], mpl.rcParams["ytick.labelsize"] = 25,25
def percent_change(series):
# Collect all *but* the last value of this window, then the final value
previous_values = series[:-1]
last_value = series[-1]
# Calculate the % difference between the last value and the mean of earlier values
percent_change = (last_value - np.mean(previous_values)) / np.mean(previous_values)
return percent_change
# Apply your custom function and plot
prices_perc = prices.rolling(20).apply(percent_change)
prices_perc.loc["2014":"2015"].plot()
plt.legend(loc="lower left", prop={'size': 25})
plt.show()
"""You've converted the data so it's easier to compare one time point to another. This is a cleaner representation of the data. """
Handling outliers
In this exercise, you'll handle outliers - data points that are so different from the rest of your data, that you treat them differently from other "normal-looking" data points. You'll use the output from the previous exercise (percent change over time) to detect the outliers. First you will write a function that replaces outlier data points with the median value from the entire time series.
Instructions:
- Define a function that takes an input series and does the following:
- Calculates the absolute value of each datapoint's distance from the series mean, then creates a boolean mask for datapoints that are three times the standard deviation from the mean.
- Use this boolean mask to replace the outliers with the median of the entire series.
- Apply this function to your data and visualize the results using the given code.
def replace_outliers(series):
# Calculate the absolute difference of each timepoint from the series mean
absolute_differences_from_mean = np.abs(series - np.mean(series))
# Calculate a mask for the differences that are > 3 standard deviations from zero
this_mask = absolute_differences_from_mean > (np.std(series) * 3)
# Replace these values with the median accross the data
series[this_mask] = np.nanmedian(series)
return series
# Apply your preprocessing function to the timeseries and plot the results
prices_perc = prices_perc.apply(replace_outliers)
prices_perc.loc["2014":"2015"].plot()
plt.legend(loc="lower left", prop={'size': 25})
plt.show()
"""Since you've converted the data to % change over time, it was easier to spot and correct the outliers."""
prices_perc
Engineering multiple rolling features at once
Now that you've practiced some simple feature engineering, let's move on to something more complex. You'll calculate a collection of features for your time series data and visualize what they look like over time. This process resembles how many other time series models operate.
Instructions:
- Define a list consisting of four features you will calculate: the minimum, maximum, mean, and standard deviation (in that order).
- Using the rolling window (prices_perc_rolling) we defined for you, calculate the features from features_to_calculate.
- Plot the results over time, along with the original time series using the given code.
prices = pd.read_csv('./datasets/prices_ebay_yahoo_nvidia.csv', index_col=0, parse_dates=True)
ebay_prices = prices[["EBAY" ]]
prices_perc = ebay_prices.rolling(20).apply(percent_change)
prices_perc
prices_perc_rolling = prices_perc.rolling(20, min_periods=5, closed='right')
# Define the features you'll calculate for each window
features_to_calculate = [np.min, np.max, np.mean, np.std]
# Calculate these features for your rolling window object
features = prices_perc_rolling.aggregate(features_to_calculate)
# Plot the results
ax = features.loc[:"2011-01"].plot()
prices_perc.loc[:"2011-01"].plot(ax=ax, color='k', alpha=.2, lw=3)
ax.legend(loc=(1.01, .6),prop={'size': 25})
plt.show()
Percentiles and partial functions
In this exercise, you'll practice how to pre-choose arguments of a function so that you can pre-configure how it runs. You'll use this to calculate several percentiles of your data using the same percentile() function in numpy.
Instructions:
- Import
partialfromfunctools. - Use the
partial()function to create several feature generators that calculate percentiles of your data using a list comprehension. - Using the rolling window (
prices_perc_rolling) we defined for you, calculate the quantiles usingpercentile_functions. - Visualize the results using the code given to you.
from functools import partial
percentiles = [1, 10, 25, 50, 75, 90, 99]
# Use a list comprehension to create a partial function for each quantile
percentile_functions = [partial(np.percentile, q=percentile) for percentile in percentiles]
# Calculate each of these quantiles on the data using a rolling window
prices_perc_rolling = prices_perc.rolling(20, min_periods=5, closed='right')
features_percentiles = prices_perc_rolling.aggregate(percentile_functions)
# Plot a subset of the result
ax = features_percentiles.loc[:"2011-01"].plot(cmap=plt.cm.viridis)
ax.legend(percentiles, loc=(1.01, .5),prop={'size': 25})
plt.show()
Using "date" information
It's easy to think of timestamps as pure numbers, but don't forget they generally correspond to things that happen in the real world. That means there's often extra information encoded in the data such as "is it a weekday?" or "is it a holiday?". This information is often useful in predicting timeseries data.
In this exercise, you'll extract these date/time based features. A single time series has been loaded in a variable called prices.
Instructions:
- Calculate the day of the week, week number in a year, and month number in a year.
- Add each one as a column to the
prices_percDataFrame, under the namesday_of_week,week_of_yearandmonth_of_year, respectively.
ebay_prices = all_prices['EBAY']
prices_perc = ebay_prices.rolling(20).apply(percent_change)
prices_perc.loc['2014-01-02':'2016-1-31']
Creating time-shifted features
Special case:Auto-regressive models
Visualize regression coefficients
Auto-regression with a smoother time series
Cross-validation with shuffling
Cross-validation without shuffling
Time-based cross-validation
Stationarity
Bootstrapping a confidence interval
Calculating variability in model coefficients
Visualizing model score variability over time
Accounting for non-stationarity